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Abstract. Differences between observed and theoretical 
eigenfrequencies of the Sun have characteristics which 
identify them as arising predominantly from properties 
of the oscillations in the vicinity of the solar surface: in 
the super-adiabatic, convective boundary layer and above. 
These frequency differences may therefore provide useful 
information about the structure of these regions, precisely 
where the theory of solar structure is most uncertain. 

In the present work we use numerical simulations of the 
outer part of the Sun to quantify the influence of turbulent 
convection on solar oscillation frequencies. Separating the 
influence into effects on the mean model and effects on the 
physics of the modes, we find that the main model effects 
are due to the turbulent pressure that provides additional 
support against gravity, and thermal differences between 
average 3-D models and 1-D models. Surfaces of constant 
pressure in the visible photosphere are elevated by about 
150 km, relative to a standard envelope model. 

As a result, the turning points of high-frequency modes 
are raised, while those of the low-frequency modes remain 
essentially unaffected. The corresponding gradual lowering 
of the mode frequencies accounts for most of the frequency 
difference between observations and standard solar mod- 



els. Additional effects are expected to come primarily from 
changes m the physics ol the modes, m particular from the 
modulation of the turbulent pressure by the oscillations. 

Key words: Sun: interior - Sun: oscillations 



1. Introduction 

In standard solar models, the stratification of the convec- 
tion zone is determined by mixing- length theory (MLT), 



hydrodynamics to a one-parameter family of models. MLT 
solar models suffer from several basic inconsistencies. For 
example, they predict that convective velocities, of several 
kms -1 , disappear abruptly in a few tens of km, imme- 
diately below the solar surface. In contrast, observations 
show convective cells with those same characteristic veloc- 
ities, and with horizontal sizes of 1000 - 5000 km, whose 
velocity fields obviously cannot vanish so abruptly. Indi- 
rect evidence from spectral line broadening indeed shows 
that the photosphere is pervaded by a velocity field with 
rms Mach numbers of the order of 0.3, and yet, in standard 
solar models these layers are assumed to be entirely quies- 
cent. Clearly, such a discrepancy between the theoretical 
description and the observations should be regarded as a 
warning not to take the quantitative predictions of the 
theory too seriously. 

Helioseismology provides quantitative diagnostics that 
pertain precisely to these critical surface layers, since this 
is where the upper turning points of the majority of modes 
are located. Thus, analysis of the observed properties of 
these modes may help clarify the consequences of the in- 
consistencies inherent in MLT and its more recent sib- 



lings (Canuto & Mazzitelli, 1991; Canuto & Mazzitelli 



1992; Canuto et al., 1996). Indeed, as we discuss in more 



detail below, adiabatic oscillations of MLT models show 
significant systematic discrepancies when compared with 
measured solar frequencies. 

It is natural, therefore, to seek to use helioseismology 
applied to these differences to improve the theoretical de- 
scription. However, this procedure is undermined by our 
present uncertainty about the physics of the oscillations 
near the top of the convection zone where they are likely to 
be strongly coupled to both the convective and radiative 
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Fig. 1. Measured frequency residuals in the sense (ob- 
servations) - (model), scaled by Q n i. The computed fre- 



quencies are for Model S of Christcnscn-Dalsgaard et al 
(1996)| The symbols indicate the source of the observed 
frequencies: + are LOWL data (Tomczyk et al., 1995) and 



o are data from Bachmann et al. (1995) 



One approach to resolving this problem has been 
the time- dependent non -local non-adiabatic mixing-length 
theory of Gough (1977) and Balmforth (1992abc) in which 
the coupling of the oscillations to both the convection and 
the ra diation is included within the framework of mixing- 



rium and is based on a simple analytical T-r relation, 
obtained from a fit to the Harvard-Smithsonian Reference 
Atmosphere (Gingerich et al., 1971). Evolution starts from 
a chemically homogeneous zero-age main-sequence model, 
and the model is calibrated to have the current solar ra- 
dius and luminosity at an age of 4.6 10 9 years. The depth 
of the convection zone in the model of the present Sun 
is O.2885i?0, close to the helioseismically inferred values 
(O.287i? ±O.OO3i?Q,~ 



Christcnsen-Dalsgaard et al., 1991 



Basu fc Antia, 1997) 



length itheory. Another approach to the problem has beeii 



O.287i? ±0.001i? Q , 

In Fig. |l| we show the differences between adiabatic 
eigenfrequencies of the standard model and meas ured solar 
p-mode frequ enci es compiled from the da ta of Tomczyk 
et al. (1995) and Bachmann et al. (1995)| . The frequency 



proposed by |Zhugzhda fc Stix (1994)| who have developed 
a model of the modal effect on mode frequencies due to 
advection of the oscillations by spatially varying radial 
flows. This approach has not yet been developed to the 
stage where it can be usefully applied to realistic solar 
models with stratification and turbulent pressure. Finally, 
Riidiger et al. (1997) have used turbulence closure assump- 
tions to parameterize the propagation of acoustic distur- 
bances through a convecting medium. 

In the present work, we use an alternative technique 
for estimating the model effects, based on the results ob- 
tained from numerical simulations of turbulent convection 
in a radiating fluid. We show that p modes can be calcu- 
lated from a mean model with hydrostatic and thermody- 
namic stratification obtained by appropriately weighted 
averages of the simulation results. We proceed by making 
simplifying, and certainly very naive, assumptions about 
the modal effects, postponing their detailed study to sub- 
sequent papers. 

We begin (Section 2) with a brief discussion of the he- 
lioseismic data and the discrepancy between measured fre- 
quencies and those calculated from MLT models. We then 
discuss (Section 3) the averaging techniques needed to an- 
alyze the radial oscillations of a convecting layer, describe 
the model computations (Section 4), and investigate the 
resulting oscillation frequencies (Section 5). Finally (Sec- 
tion 6), we discuss the relevance and limitations of the 
results, and indicate future plans. 

2. Frequency Residuals for a Standard Solar 
Model 

To illustrate the differences between the observed frequen- 
cies and those of current "standard" solar models we con- 



differences have been scaled with the quantity Q n i defined 
as the ratio of the mode mass of a mode with radial order 
n and degree / to the mode mass of a radial mode of the 
same frequency, the mode mass being defined by 



\£\ 2 p dV = M mod0 |£(i? G 



(1) 



where £ is the displacement eigenfunction of the mode. 
Here we have introduced the equilibrium density stratifi- 
cation po and the integral is over the volume of the Sun. 

Perturbations concentrated in the surface region of the 
Sun give rise to frequency residuals which, when scaled 
with Q n i, are a function of frequency except at high de- 



grees (e.g. Christensen-Dalsgaard k Berthomieu, 1991); 
it is evident that this is, indeed, the dominant trend in 
Fig. |l|. We therefore conclude that the principal contribu- 
tion to the solar frequency residuals arises from the surface 
layers. This is consistent with the theoretical arguments 
that these layers are particularly badly represented in the 
models. 

Some insight into the effects of this region on the fre- 
quencies might be obtained by linearizing the relation be- 
tween the structure and the frequencies in terms of dif- 
ferences between the Sun and a reference model. To the 
extent that the oscillations can be regarded as being adi- 
abatic, the frequencies are determined by the hydrostatic 
structure of the model as well as by the adiabatic expo- 
nent Ti = (d hip/d lnp) s relating pressure p and density 
p, the derivative being at constant specific entropy s. For 
the purpose of analyzing near-surface effects a convenient 
variable is v = Ti/c = 2cj c /g, where c = (Tip/p) 1 / 2 is the 
adiabatic sound speed and tu c is the isothermal acoustic 



cutoff frequency. It was shown by Christensen-Dalsgaard 



is computed with the OPAL equation of state (Rogers 
et al, 1996| ), OPAL opacities ( |Iglesias et al. ,1992| ) ; ~ 

tP7 ( 1 QQ1 M nna r>i f \ 



sider lY fr odcl S of Christensen-Dalsgaard et al. (1996)|. This fc Thompson (1997) that if the differences are expressed 



in terms of Lagrangian differences (i.e., differences at fixed 
mass) S m v/v and S m c/c, the contribution from 5 m c/c is 
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Fig. 2. Scaled kernels relating the Lagrangian change in 
v = Ti/c to frequency changes (cf. Eq. |^), plotted against 
frequency, for modes of degree / < 5. Kernels are shown at 
the photosphere (solid curve) and at the following depths 
below the photosphere: 2 10~ 4 i? (dotted curve), 5 10 _4 i? 
(dashed curve), 10 _3 i? (dot-dashed curve), and 5 10 _3 i? 
(triple-dot-dashed curve) 



where the kernels K™ l c (r) can be calculated from the struc- 
ture and eigenfunctions of the reference model. 

Fig. U shows —Q n iK™ l c (r) for modes of degree I < 5, 
plotted against frequency at several depths; this corre- 
sponds to the effects on the frequencies of a (negative) 
delta-function change located at the selected depths. The 
scaling with Q n i was included to suppress the effect of 
the varying mode mass, whereas the sign is changed to 
ease comparison with the differences between observations 
and model in Fig. [j]. Two features of the kernels are im- 
mediately evident: they are very small at low frequency 
and they become increasingly oscillatory with increasing 
depth of the perturbation to the model. Both features 
may be understood quite simply in terms of the prop- 
erties of the oscillation eigenfunctions. At low frequencies 
the upper turning points of the modes are located at rela- 
tively greater depth; hence the eigenfunctions are evanes- 
cent in the region considered and the kernels are there- 
fore small. At higher frequency the modes are oscillatory 
quite close to the surface; as frequency is varied the nodes 
of the eigenfunction move relative to the location consid- 
ered, giving rise to the oscillation, the shift being more 
rapid at greater depth. We note that the observed differ- 
ences in Fig. |l| are indeed very small at low frequency and 



3. Averaged Equations for Radial Oscillations 

Our analysis of the equations of momentum and conti- 
nuity will closely follow that given by Stein & Nordlund 



1991)| in the context of their analysis of mode excitation. 
However, the treatment of the energy equation will require 
some additional assumptions in order to obtain a closed 
form for the resulting equations. 

We begin by introducing a horizontal (but not tempo- 
ral) average denoted by (...) and defining 

_ (PU Z ) 



(P>> Pg = (Pg). and 



(p) 



(3) 



for the mean density, gas pressure, and vertical velocity. 
We then define fluctuations around these quantities by 



p = p — p etc. 



(4) 



The quantities denoted by the overbars thus contain both 
the mean stratification and the p-mode oscillations whilst 
the primed variables are convective fluctuations. Follow- 



ing Stein & Nordlund (1991) it can be shown that the mo- 
mentum equation (neglecting dissipative terms) becomes 
exactly 

dpu 
~dt 

The quantity (pu'^) is the turbulent pressure which we 
denote by pt- The total pressure is then defined by 



°[pui + (pu?)]-^+gp 



(5) 



(6) 



(7) 



P = Pg + Pt ■ 

We now linearize Eq. (^) by writing 

p = p (z) + epi(z,t) , 

with similar expressions for p and u z with the proviso that 
u z0 = 0, i.e., that there is no net mass flux through any 
horizontal surface. Thus we obtain the linearized momen- 
tum equation 

du zl dpt 

Pa ^r^-^ +9Pl - (8) 

Similarly one may show that p and u z satisfy the usual 
continuity equation 

<9p + dpu z 
dt 



dz 



o 



(9) 



The energy equation is 

Dp g _ PgFi Dp 
Dt 



= - (r 3 - 1) v • F iad 



(10) 



p m 

where Ti and T3 are the usual thermodynamic variables 
and .Frad is the radiative heat flux. Averaging as above 
and using equation (9) yields the equation 

drTrr dv„ . . du-, 
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In order to close the system we must make a number 
of approximations. The first term on the right-hand side 
of Eq. (FyJ) represents the compressive work and has pre- 



viously been argued to be negligible (Stein & Nordlund 
199lj)7]For a perfect gas, dividing p s by T3 — 1 yields the 
energy so the second term on the right hand side is the 
convective flux divergence multiplied by (F 3 — 1) while 
the third term is evidently the radiative flux divergence 
multiplied by the same factor. Their sum is therefore pro- 
portional to the total flux divergence if departures from a 
perfect gas law can be ignored. 

In order to proceed with analytical methods to arrive 
at an equation for the pulsations, we now make the sim- 
plest possible assumption in Eq. (fill), namely that the 
time- varying part of the sum of the right-hand-side terms 
may be neglected. This is obviously just a first, simplified 
step in a more complete investigation of the interaction 
between the convective fluctuations and the oscillations; 
even so, it is of interest to investigate the properties of the 
oscillations if the fluctuations of the 3-D averages vanish 
and compare the outcome with the observed frequencies. 

We note, however, that the validity of the assumption 
is clearly suspect since the time scales for convective mo- 
tion are comparable with the periods of the oscillations. In 
the future we will address the validity of this assumption 
quantitatively, using direct measurements of the ampli- 
tudes and phases of these terms in numerical simulations. 

In relating the variation in p g to the variation in total 
pressure the two simplest assumptions are 

1. that the Lagrangian variation in the turbulent pressure 
may be neglected, or 

2. that the Lagrangian variation in the turbulent pressure 
is directly proportional to the Lagrangian variation of 
the gas pressure. 



In the first case, 



Di 



Dp 

dT 



(12) 



leading to the linearized energy equation 

dpi . dp - r du zl 
-^r +u zl — + Ttfo—— = , 
at oz oz 

where we have introduced the so-called "reduced gamma" 
defined by 



(13) 



(?>g r i)o 



Po 



(14) 



In the second case, if the relative variation is the same in 
both components of the pressure, the effective gamma is 
the average value (Li)o for the gas. 

The reduced gamma was originally introduced by 



In the present work we investigate the influence of the 
response of the turbulent pressure by using two somewhat 
extreme cases: the reduced Li (Eq. |l4]) and average gas 
Li. The choice of Li thus constitutes our description of 
the modal effects contributing to the measured frequency 
residuals. The effects of model discrepancies are dealt with 
by using a numerical simulation of the surface layers to 
replace the usual MLT prescription, as we discuss in detail 
in the next section. 

Equations (ph, (0) and ([l3|), when supplemented with 
a prescription for Ti, now constitute a closed set of equa- 
tions for the radial oscillations of the convecting layer. 
Their form is identical to the equations for linear, radial 
oscillations of a hydrostatic layer defined by a pressure 
variation po(z) and an adiabatic exponent Ti given by ei- 
ther f \(z) or (ri)o. In order to use results of a simulation 
to determine the lowest-order effect on the frequencies of 
oscillations it is therefore only necessary to specify these 
two mean quantities. 

We finally note that, even though the derivation pre- 
sented here is formally only valid for radial oscillations, 
we shall apply the formalism to nonradial oscillations as 
well; as in the radial case, this is done by replacing Ti 
in the oscillation equations by either T\(z) or (Li)o. For 
modes of relatively low degree, where the motion is pre- 
dominantly vertical in the near-surface layer, this is likely 
to be a good approximation, while it is more questionable 
for high-degree modes where the horizontal and vertical 
components of displacement are comparable. 

4. Model computations 

Our principal goal is to investigate how the treatment 
of convection affects the computed structure and oscil- 
lation frequencies in models that cover a large fraction of 
the Sun. To achieve this we consider averages of hydro- 
dynamical simulations, extended continuously with enve- 
lope models computed using the MLT. We compare these 
patched models with purely MLT envelopes using other- 
wise the same physics, with the observed frequencies, as 
well as with the properties of standard evolutionary solar 
models such as Model S considered in Section 2. 

Detailed descriptions of the simulations have been pro- 
vided by Ste in & Nordlund ( |1989| ; |1998| ) and [Nordlunc 1 
|st al. (1996) . It is a 3-dimensional calculation, incorpo- 
rating LTE radiative transfer and provisions for using an 
arbitrary, tabular equation of state. Most of our results are 
based on a simulation corresponding to a rectangular box 
of horizontal dimensions 6 x 6Mm and a vertical extent of 
3.4 Mm, of which 2.9 Mm are in the convectively unsta- 
ble region. This was resolved with a mesh with 100 x 100 
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Fig. 3. A vertical slice through the simulation showing the temperature variation. Also shown are the horizontally 
averaged superadiabatic gradient and ratio of turbulent to gas pressure 



et al., 1988; Mihalas ct al., 1990| ), and an updated version 
of the Kurucz (1991) opacities (cf. Trampcdach, 1997). 
The chemical composition corresponded to abundances by 
mass of hydrogen and heavy elements of X = 0.736945 and 
Z = 0.0180550, respectively; this is close to the envelope 
composition of Model S. The model assumed a gravita- 
tional acceleration of 2.740 10 4 cms~ 2 and an average ef- 
fective temperature of 5777 K. The averaged models were 
constructed by averaging first in the horizontal direction 
and subsequently carrying out a Lagrangian time average, 
i.e., a time average on a fixed mass scale, to filter out the 
main effect of the p modes excited in the simulation. The 
models were averaged over at least half an hour of solar 
time. Figure || shows a vertical slice through a single time 
step of the calculation. The plot also indicates the aver- 
aged superadiabatic temperature gradient and the ratio 
of turbulent pressure to gas pressure. Evidently there is a 
narrow superadiabatic layer near the top of the box and 
the turbulent pressure is significant only near this region. 
At greater depths in the box the stratification is nearly 
adiabatic. 

To test the influence of numerical resolution, as well 
as of the handling of radiative transfer etc., we have car- 
ried out a number of, generally shorter, simulations with 
varying numbers of mesh points (experiments drawn upon 
in the present work have resolutions 32 2 x 41, 50 2 x 82, 
63 2 x 63, 100 2 x 82, 125 2 x 82, 125 2 x 163, and 253 2 x 163— 
some of these were calculated with earlier versions of the 
equation of state, and have slightly different chemical com- 
positions). 

Since the simulation covers only a narrow region near 
the top of the solar convection zone, we require an enve- 
lope to which it can be matched. This was constructed 
using an MLT envelope code. It used the same (MHD) 
equation of state as did the simulation; also, the hydro- 
static equilibrium in the envelope calculation included a 
turbulent pressure of the form 



Pt = flpv, 



2 

con 



(16) 



The simulation and envelope models were patched to- 
gether to produce smooth 1-dimensional models. This was 
achieved by adjusting (3 and the mixing length in the MLT 
treatment such that, at a fixed pressure near the bottom 
of the simulated region, turbulent pressure, temperature 
and density matched continuously. The kinetic energy flux 
was neglected in the envelope part of the model; its inclu- 
sion would change the superadiabaticity of the envelope 

oln nrl-iFl i r V-vt n F F Vi a £ifl-<a/"tF nrAi 1 1 /-I V-\/a /-^ m F d cm nil Din/ifl FV\ci an 



adiabatic part of the convection zone corresponding to 
the simulation. In particular, it is possible to determine 
the depth of the convection zone predicted by the adiabat 
obtained in the simulation. Strikingly, the 100 2 x 82 simu- 
lation considered here leads to a convection-zone depth of 
O.286i?0, very close to the helioseismically i nferred value 
(0.287i? Q ± 0.001i? Q , |Basu fc Antia, 1997| ). In view of 



the fact that no parameters were adjusted to achieve this 
agreement, this constitutes additional evidence of the con- 
sistency of the model. 

It should be noticed that the envelope model does not 
allow for convective "undershoot" below the convection 
zone proper. Thus, either the undershoot is indeed small, 
or the agreement of the convection zone depth is fortu- 
itous. There is independent evidence, from simulations of 
a convection zone with greatly enhanced total flux (Stein 
et al. 1997, Dorch 1998) that the convective undershoot is 
actually small. 

At a total flux ~ 3 10 5 solar, the convective undershoot 
is of the order of 10 - 20 % of the solar radius. The under- 
shoot may be expected to scale in rough proportion to the 
velocity (Hookes law), which again scales roughly as the 
third root of the total flux. The undershoot at the nominal 
solar flux is therefore expected to be only a fraction of a 
percent of the solar radius. 

Errors in the equation of state, or uncertainties in the 
precise abundances of chemical elements could also affect 
the apparent match between the predicted and observed 
depths of the convection zone. 

4-1. Computation of p-mode frequencies 

Two types of patched models were considered when com- 
puting frequencies of adiabatic oscillations. In one, in the 
following referred to as reduced-gamma models (RGM), 
the adiabatic exponent between pressure and density was 
replaced by the reduced Ti defined by Eq. (14); in the 



second, the gas-gamma models (GGM), the adiabatic ex- 
ponent was obtained as the average (ri)o- This distinction 
affects the determination of the adiabatic sound speed c, 
v and evidently the oscillation frequencies. These choices 
of Ti are certainly naive, and we do not expect to obtain 
a perfect match to the observed frequencies with any one 
of them. However, as stated in the Introduction, in the 
present paper we are mainly concerned with the model ef- 
fects; i.e., we wish to explore the effects on the oscillation 
frequencies that come purely from changes in the mean 

c+vaf ifir*af inn 
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Fig. 4. Logarithmic Lagrangian (at fixed mass) differ- 
ences in v = Ti/c between patched and comparison en- 
velope models (solid and dashed lines) and between the 
GGM envelope and standard solar model (Model S) 



with the Rosseland opacity corresponding to the non-grey 
opacity used in the simulations. Furthermore, the atmo- 
spheric structure was based on a T-r relation obtained 
as an analytical f it to the mean T-r relatio n for the sim- 
ulation (see also Trampedach et al., f999 ). The mixing- 
length parameter was adjusted to give exactly the same 
convection-zone depth as in the corresponding matched 
models. 

Figure [| shows the difference in the quantity v = 
yjTipo/po between the GGM and SEM and between the 
RGM and SEM. The differences were evaluated at fixed 
mass (or, equivalently, at fixed hydrostatic pressure). In 
the case of the GGM, therefore, the difference largely re- 
flects the difference in p , caused by the presence of tur- 
bulent pressure and the consequent reduction in the gas 
pressure; there is an additional, comparatively small, con- 
tribution arising from the change in Ti in the region of 
steeply increasing hydrogen ionization. The combined ef- 
fect is a reduction of v by up to about 15% in a narrow 
region centered on the superadiabatic layer. In the RGM, 
v is further reduced by the replacement of (ri)o by the 
reduced f J. Note that the interior part of the two matched 
models and the SEM were calculated with slightly differ- 
ent mixing-length parameters (as the SEM model has no 
turbulent pressure whereas the matched envelope models 
include a turbulent pressure according to Eq. (fl6|)), so the 
models are not exactly identical below the matching point. 
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Fig. 5. Scaled differences between frequencies of patched 
and comparison envelope models (crosses) and between 
the GGM envelope and standard solar model (Model S) 
(diamonds) 



5. Results on oscillation frequencies 

5.1. Comparison of frequencies 

Adiabatic oscillation frequencies were calculated for the 
patched and comparison envelopes, using the procedures 
described by Christcnsen-Dalsgaard fc Bcrthomicu (1991)| . 



The models are not complete evolved solar models, unlike 
the standard solar model used in the construction of Fig. |[ 
Nevertheless we believe that the differences in mode fre- 
quencies between the models will be an accurate represen- 
tation of the actual change in the mode frequencies of a 
standard solar model constructed under the same physical 
assumptions. In particular, by restricting the comparison 

f A ^Vli~vQ£l TYIA^DC TTrtlir>Vl OTO I 1~ > I 1 \ T \ Z ] in + Vl C\ ^rtTlirfyif 1 An '//111 £\ 
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Fig. 6. Measured frequency residuals in the sense (obser- 
vations) - (model), scaled by Q„;. The model frequencies 
are for the GGM, with the gas IV The symbols indicate 
the source of the observed frequencies: + are LOWL data 



1.000 



r/R 
0.999 



0.998 



et al. (1995) 



(Tomczyk et al., 1995) and o are data from Bachmann 



inclusion of settling the structure of the radiative interior 
of the latter model differs substantially from the structure 
of the envelopes. 

Scaled frequency differences between the patched and 
comparison envelopes are shown in Fig. ||. To avoid clut- 
tering the diagram only modes with I < 300 were included. 
For these, the scaled differences between the envelope fre- 
quencies are indeed predominantly functions of frequency 
which become very small at low frequency. Also, interest- 
ingly, the difference between the GGM and SEM frequen- 
cies is of a similar magnitude and shape as the difference 
between the observed and Model S frequencies shown in 
Fig. [I], the GGM frequencies being decreased by up to 
about 15/iHz relative to the SEM frequencies. On the 
other hand, in accordance with Fig. ^, the frequency dif- 
ferences for the RGM frequencies are substantially larger. 

To interpret Fig. |l], a more appropriate comparison is 
evidently between the frequencies of Model S and those 
of the models including the averaged simulations. Accord- 
ingly, Fig. H also shows scaled differences between the fre- 
quencies of the GGM and Model S. These are slightly 
smaller than the corresponding differences between the 
GGM and SEM; also, they show somewhat larger scatter, 
as a result of model differences extending over a larger 
part of the convection zone. 

As argued above, the near agreement between the 
depth of the convection zone in the patched models and in 
the Sun indicates that the structure of the deeper parts of 
the convection zone may be in similar agreement. Thus it 
is of evident interest to compare the computed frequencies 
with the observed values. A comparison of Figs and | 
strongly indicates that the GGM is closer to the solar data 
than is the RGM. Accordingly, in Fig. || we show scaled dif- 
ferences between observed frequencies and those of GGM, 
again restricted to modes with vj(l + 0.5) < 50 ^Hz and 
hence trapped in the convection zone. The general magni- 
tude of the p-mode frequency differences is substantially 
reduced, and the differences show even less dependence 
on degree than in Fig. |], indicating that the origin of 
the remaining differences is concentrated very close to the 
surface. Note also that differences in Fig. |l| decrease sig- 
nificantly from zero at rather higher frequency than do 
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Fig. 7. Pressure as a function of depth for an averaged 
3-D model (full drawn), and for the comparison standard 
envelope (SEM) model (dashed). The dashed-dotted curve 
shows the pressure stratification of a 3-D model where 
the gradient of the turbulent pressure has been artificially 
removed from the vertical pressure balance. The upper 
abscissa shows the corresponding position in a complete 
model, in terms of the fractional radius r/R 



to the surface than are the differences between GGM and 
Model S. 

The gas Ti simulation is rather successful in repro- 
ducing the observed mode frequencies. In contrast, the 
reduced Ti gives rise to a frequency shift which is substan- 
tially higher than observed. Mode physics effects probably 
account for the remaining discrepancies in Fig. |(| These 
change sign as a function of frequency, indicating that the 
effects producing them depend on depth or frequency or 
both. The frequency differences for the f modes, which 
are barely affected by changes in the hydrostatic struc- 
ture of the model or in Fx, are essentially the same as 
in Fig. and of similar magnitude to the other modes 
now. Thus, there must be additional mode physics effects 
beyond those which can be included in a Ti(z,u>). 

5.2. Effects of turbulent pressure 

The turbulent pressure pt = (pu®) enters the momentum 
equation (||) in an exactly analogous manner to the gas 
pressure. It thus adds directly to the gas pressure and 
provides additional support against gravity, resulting in 
an elevation of the solar surface. The extent of the el- 
evation can be determined by comparing the heights of 
comparable surfaces of pressure, density or acoustic cutoff 
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envelope (SEM) model. This model has, by definition, the 
same average T(r) relation as the 3-D model, and uses the 
same opacities, and a conventional mixing length recipe 
with the mixing length chosen such as to give the same 
asymptotic adiabat at depth. 

The total elevation of the 3-D model relative to the 
standard envelope model is about 150 km in the mid 
photosphere, and sets in already slightly below the pho- 
tosphere. The average temperature in the 3-D model is 
higher than that of the 1-D model in the surface re- 
gions, and the correspondingly higher pressure scale height 
causes an elevation that adds to the one caused by the tur- 
bulent pressure. 

The difference in thermal structure reflects fundamen- 
tal differences between 3-D and 1-D models, rather than 
just differences in the efficiency of convective energy trans- 
port; a 1-D model with a temperature structure consistent 
with the average pressure stratification of the 3-D model 
would have too high an effective temperature, as mea- 
sured by the emitted surface radiation, while the radiation 
emitted from the 3-D model is consistent with the solar 
effective temperature. The difference is caused by the tem- 
perature sensitivity of the opacity; relatively hot regions 
are more opaque, and hence contribute less to the emitted 
radiation, while relatively cold regions are less opaque and 
hence contribute more. The higher mean temperature of 
the 3-D model is thus "hidden" from view, but is reflected 
in the pressure stratification. 

The total photospheric elevation illustrated in Fig. [7] 
provides a simple if crude understanding of the frequency 
changes found for the GGM model. The high-frequency 
modes have upper turning points in the photosphere; thus 
the location of the upper turning point is lifted by the 
elevation and hence the acoustic cavity is extended, lead- 
ing to a reduction in the frequency. For radial modes the 
relative frequency change resulting from this effect can be 
estimated as 

Sv _ Ar/c ph 

where Ar is the elevation, c p h is the photospheric sound 
speed and the denominator is the acoustical radius tq 
of the star. Adopting an elevation Ar = 150 km, c p h = 
7kms _1 and using To ~ 3700s, we find 8v/v ~ 610~ 3 , in 
general agreement with the frequency changes shown in 
Fig-!- 

This argument is roughly consistent with the expres- 
sion for the frequency change given in Eq. (|^). Using the 
asymptotic properties of the cigenfunctions, neglecting the 
details of the behaviour near the upper turning point, Eq. 
(§) may be approximated by 




P 9 [cgs] 

Fig. 8. Mean turbulent pressure divided by gas pressure 
as a function of depth for three different runs, with dif- 
ferent numerical resolutions: 253 2 x 163, 125 2 x 163, and 
63 2 x 63. 



for the right-hand side of Eq. ( [l8| ) 

f R S m v dr ^ I f R 5 m p dr ^ 1 x f R <W dr 
Jo vc~2j p c ~ 2 C P h J p 

- l<£8 m r , (19) 

where we took c ~ c p h to be constant over the region 
affected by turbulent pressure and, in the last equality, 
used the equation of mass conservation. Identifying the 
Lagrangian radius change S m r with Ar we recover Eq. 
(|l7|), apart from the factor 1/2. 

5.3. Effects of numerical treatment 

The ratio of turbulent pressure to gas pressure is shown 
in Fig. ||. The figure shows that the shape of the relation 
is a very robust property, but that there is an overall scale 
factor that depends slightly on the numerical resolution. 
The question thus arises whether the turbulent pressure 
has converged and whether it is possible to obtain inde- 
pendent, observational constraints on its convergence. 

The turbulent pressure peaks in a layer that is essen- 
tially at the solar surface, and the velocity distribution 
that contributes to the turbulent pressure also produces 
the observed (macro- and micro-) turbulent broadening of 
photospheric spectral lines. Figure || shows that the lo- 
cation, shape and width of the maximum in the ratio of 
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Fig. 9. Comparison of observed and synthetic Fell A5414 
spectral line profiles, based on the time average of syn- 
thetic spectral line profiles from a 253x253x163 simu- 
lation. The observed spectral line (from the Liege atlas, 
Dclbouille et al., 1973| ) is shown full drawn, while the syn- 
thetic spectral line is shown with (diamond) symbols. A 
synthetic spectral line calculated with u z = is also shown 
(dashed line with + symbols). 



The width of weak iron lines observed at solar disc cen- 
ter is a direct measure of the vertical velocity amplitudes, 
and is thus an excellent proxy for the turbulent pressure. 

Iron lines are ideal diagnostics of photospheric dynam- 
ics, not only because iron is a sufficiently heavy element 
for the thermal Doppler speed to be smaller than typical 
photospheric velocities, but also because of the large num- 
ber of iron lines that have accurately measured absorption 
coefficients and wavelengths. Many photospheric lines are 
blended, but because of the large number of iron lines, it 
is still possible to find dozens of relatively clean lines of 
Fel and Fell in the solar spectrum. 

Weak Fell lines have the advantage that they are 



the velocity; we thus estimate that the turbulent pressure 
magnitude of the model is within about 2% of the correct 
value. 

The similarity of the density stratifications for the var- 
ious runs corroborates this argument. Even between our 
lowest and highest resolution the change in the elevation 
would hardly be visible in Fig. [?]. Hence, we expect that 
the elevation of the atmosphere is nearly fully accounted 
for at our highest resolutions. 

In addition to the turbulent pressure stratification, the 
thermal structure is also important in determining the 
mode frequencies. The thermal stratification itself is quite 
insensitive to the numerical resolution — there is not much 
difference even between a model with a resolution of only 
32x32x41 and one with a resolution of 253 x 253 x 163 (cf. 
[Stein fc Nordlund, 199S| , Fig. 26). The robustness of the 
average temperature structure is partly due to the strong 
constraints imposed by stratification (Stein & Nordlund 
1989, Spruit et al. 1990, Spruit 1997, see also Stein & 
Nordlund 1998, Fig. 13). 

"Spectral line blocking" plays a critical role in deter- 
mining the surface temperature (cf. Mihalas, 1978| , pp. 
167-169). The blocking of photons by a large number of 
spectral lines in the solar spectrum forces an increase of 
the continuum temperature by about 3%, relative to a 
case with no line blocking, in order to maintain the solar 
effective temperature. The surface temperature in turn, 
through the extreme temperature sensitivity of the opac- 
ity, has a strong influence on the surface pressure. The 
temperature and pressure determine the surface entropy, 
and are hence of direct importance for the entropy of the 
bulk of the convection zone. 

We have made a substantial effort to ensure that the 
spectral line blocking is handled as accurately as possible, 
given the approximations necessitated by the constraints 
on computer capacity. During the simulations, the spec- 



tral lines are represented by only four bins (cf. Nordlund 



formed quite close to the layer where_the ratio of tur- 1982; Nordlund & Dravins, 1990; Trampedach, 1997), but 



bulent to gas pressure peaks (cf. Fig. |8J), while weak Fel 
lines are formed about half a pressure scale height further 
up. 

Figure || illustrates the close match that is obtained 
between synthetic and observed weak Fell lines, and also 
shows what a weak Fell line would look like if there were 
no macroscopic photospheric velocities. The line is much 
to deep and narrow, and it is perfectly symmetric, in con- 
trast to the slightly asymmetric observed spectral line. On 
the other hand, when the full velocity and temperature 
fluctuations of the 3-D models are included, and a syn- 
thetic spectral line is computed as an average over space 
and time, the resulting line matches the width and shape 



in n ncn - 



when the lookup tables for the binned opacities and source 
functions are computed we ensure, by calibrating against 
a detailed, 2748 frequency point calculation for a 2-D slice 
from the model, that the resulting surface flux is correct, 
i.e., that the net effect of the line blocking is accurately 
accounted for. 

The accuracy of the spectral line blocking is illustrated 
by the close agreement obtained between the predicted 
depth of the solar convection zone and the depth deter- 
mined fro m solar observat ion. Less complete spectral line 
data sets ( Bell et al., 1976 ) that were used in earlier work 
resulted in a spectral line blocking that was about 40% 
smaller. As a result, the surface temperature was about 
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pressure (and hence entropy) stratification of the surface 
layers; the convective efficiency of the superadiabatic lay- 
ers immediately below the surface that determines the 
transition to the asymptotic adiabat; and the equation 
of state that determines the further run of temperature 
versus pressure through the convection zone. 

In this connection it should be noted that the treat- 
ment of viscous and diffusive effects (Stein & Nordlund 
1998, Eqs. 4-5) is a matter of relatively minor importance. 
Test runs have shown that reasonable changes to the nu- 
merical coefficients have only minor effects (too small to 
be visible in figures such as Fig. 26 of Stein & Nordlund 
1998), and also do not cause visible differences in the syn- 
thetic spectral line profiles (Fig. ^). In effect, such changes 
correspond to slight changes in the numerical resolution. 
In tests with driven, isotropic 3-D turbulence, this type of 
numerical viscosity is seen to lead to an inertial range that 
extends to within about a factor of five from the Nyquist 
wave number, with an average viscous dissipation that ap- 
proaches a limiting value with increasing numerical reso- 
lution. Since, furthermore, the average viscous dissipation 
is small compared to other terms in the energy equation, 
the direct influence of uncertainties in viscous dissipation 
is not of great concern. An accurate treatment of the ra- 
diative energy transfer, on the other hand, is of central 
importance for the energy balance of the crucial surface 
layers. 

Not surprisingly, we conclude that the main factor of 
importance in the numerical treatment is the numerical 
resolution; the vertical resolution must be sufficient to 
resolve the sharp drop in temperature at the surface of 
granules, and the horizontal resolution must be commen- 
surable, to avoid very skew cell aspect ratios. Ultimately, 
whether the numerical resolution is "sufficient" or not 
must be judged by empirically studying the convergence 
of the solution, and by the consistency of diagnostics such 
as the synthetic spectral line profiles shown in Fig. ||, the 
depth of the solar convection zone, etc.. 

Such a procedure can only show consistency, and does 
not prove that the model is converging towards a correct 
solution. It is possible, for example, that the numerical 
models still contain inaccuracies that do not significantly 
affect the depth of the convection zone, but which could 
still have a noticeable effect on the surface structure. How- 
ever, the more diagnostics that are found to be consistent 
with observations, the smaller is the margin for remain- 
ing inconsistencies. The spectral line width data (Fig. ^), 
for example, does not leave much room for inaccuracies 
in the velocity amplitudes. Also, since the convective flux 
depends on the product of velocity and temperature fluc- 
tuations, errors in the magnitude of the temperature flue- 



Details could still be improved, though. A more de- 
tailed binning of the opacity could, for example, lead to a 
slightly different steepness of the sharp temperature drop 
at the surface, while leaving the overall cooling of gas 
visiting the surface almost unchanged. It will be neces- 
sary carefully to investigate this and similar issues when 
studying the more subtle mode physics effects, but for the 
purpose of the present study we are content with the nu- 
merical accuracy - the overall magnitude of the elevation 
caused by the combination of the turbulent pressure and 
the thermal differences is a result that does not depend 
on such subtleties. 

5.4. Comparison with evolution model 

We have concentrated on the direct effects of convection 
on the model and frequencies, by comparing the averaged 
hydrodynamical models (GGM and RGM) with the SEM 
model computed with mixing-length theory but otherwise 
using as far as possible the same physics as the hydro- 
dynamical models. However, it is of evident interest to 
compare also with Model S, which can be taken as repre- 
sentative of 'standard' evolutionary models of the present 
Sun. As seen from Fig. ^ by comparing the solid and dot- 
ted lines, there are in fact considerable differences between 
the SEM and Model S. We have identified the dominant 
causes of these differences as arising from two aspects 
of the treatment of the atmosphere of the models: the 
assumed T(t) relation and the low-temperature opacity. 
Both change the temperature structure, as a function of 
mass, in the atmosphere and uppermost part of the con- 
vection zone; this in turn affects the density and, in the re- 
gion of partial hydrogen ionization, Ti, leads to the model 
differences illustrated (see also Trampedach et al. 1999). 
For modes with frequencies below about 3000 /iHz the up- 
per turning point is essentially below the region where 
these model differences are significant; hence, in Fig. ^, 
the frequency differences for GGM - Model S and GGM 
- SEM are very similar. At higher frequency, however, the 
increased model differences very near the surface, when 
the SEM model is used as reference, lead to a significant 
increase in the magnitude of the frequency differences. 

6. Discussion 

In the classical theory of solar structure, a one-parameter 
family of models (MLT) is calibrated against the known 
radius of the Sun. It is well known that MLT is based on 
a number of fundamentally inapplicable and inconsistent 
assumptions, and so a large number of alternative models 
of stellar convection have been proposed. 
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between upflows and downflows as found in the simula- 
tions. The models of Lydon et al. ( 1993b] , 1993a) are 



essentially attempts to parameterize a wide range of con 
vection simulations using a formalism similar to MLT but 
incorporating also the contribution of the kinetic energy 
flux to the flux-balance equation. Canuto & Mazzitclli 
(1991)|have taken an approach based on modern theories 
of turbulence, and have attempted to produce a parame- 
terized expression based on such theories. Finally, Canuto 
( |1992| ; |1993| ) has produced an ambitious model of convec- 
tion based on a Reynolds' stress formalism. 

These more elaborate models of convection are po- 
tentially very useful, if it can be shown that they cap- 
ture essential aspects of the full 3-D convection in terms 
of much simpler equations. In particular, non-local and 
time-dependent models of convection such as the ones by 



Gough (1977), Buchler (1993), and Houdek (1997) would 



be very useful, for example in the modeling of pulsating 
stars, since full 3-D simulations are too expensive to be 
used in that context. The most obvious way to proceed 
is by use of numerical simulations such as the one used 
here or those of Kim et al. (1995), treating the simula- 



tions as data against which the models are to be tested 
and calibrated. 

A second approach to testing the simplified models, 
and one which has been more widely adopted so far, is 
the helioseismic approach in which the frequencies of os- 
cillations of a solar model constructed with a given theory 
are compared with the observed frequencies. Such a proce- 
dure, however, must deal with the difficulty that the com- 
plete structure of the surface layers is certainly underde- 
termined by the (low- and intermediate-degree) oscillation 
frequencies since, as Christensen-Dalsgaard & Thompson 
(1997)| have shown, the near-surface contribution to the 
frequencies depends only on v and not separately on, e.g., 
Pq and Ti. Thus improved agreement with measured mode 
frequencies cannot, by itself, be taken as evidence that a 
given model of convection is a better description of reality. 

In the present work we have attempted to improve on 
this approach by analyzing the problem in such a way as 
to obtain a physically justifiable description of the oscilla- 
tions. Moreover, by using a numerical simulation of con- 
vection we give ourselves no free parameters with which 
to calibrate our model. Given this approach, the fact that 
we are able to construct a complete solar envelope model 
with essentially the correct convection-zone depth is, in 
itself, a considerable achievement for the simulation. 

We have here investigated model effects on the mode 
frequencies, primarily those due to changes in pressure 
support of the atmosphere and 3-D radiation transfer. 



an inevitable consequence of the temperature dependence 
of the opacity. 

One might attempt to estimate the elevation effect 
from turbulent pressure using a local convection model 
such as MLT or the model of Can uto fc Mazzitclli (1991 , 
1992). However, as discussed by Antia fc Basu (1997) 



such models cannot accurately account for the effect of 
turbulent pressure because they result in an artificially 
abrupt upper boundary of the convection zone and there- 
fore a serious overestimate of the turbulent-pressure gra- 
dient there. In addition, as discussed in section 5.2, a 1-D 



model with the correct pressure stratification unavoidably 
has a surface radiation flux that corresponds to an incor- 
rect effective temperature. Thus, there are inherent limi- 
tations in simplified 1-D models of convection. 

The principal remaining uncertainty in determining 
the oscillation frequencies from 3-D models lies in the un- 
certain mode physics. In particular, non-adiabatic effects, 
and the response of the turbulent pressure to the com- 
pression and rarefaction in the oscillations needs to be 
understood. We have attempted to quantify this uncer- 
tainty by calculating two models: the RGM in which the 
turbulent pressure is assumed to be unaffected by the os- 
cillations and the GGM in which it is assumed to respond 
in exactly the same way as does the gas pressure. The 
frequency discrepancies for the RGM are almost exactly 
twice those for the GGM. The GGM produces frequencies 
that are closer to those observed, but this should of course 
not be taken as evidence that the turbulent pressure re- 
sponds in exactly the same way as the gas pressure. 

The actual depth- and time-dependent mode re- 
sponse of the turbulent pressure produces, together with 
the response of the gas pressure, a complex-valued, 
and frequency-dependent ri(z,o;). At any particular fre- 
quency, the real part of Ti may be expected to have a dif- 
ferent depth dependence than that assumed in both the 
RGM and GGM models; therein lies an essential part of 
the modal effects, and thus a potential explanation for 
part of the remaining differences between the observed 
and calculated oscillation frequencies. 

Preliminary investigations of the relation between 
<51n(p) and i51n(P) in numerical 3-D simulations overlaid 
with initial radial eigenmodes show that non-adiabatic ef- 
fects are indeed significant. The effective gamma appears 
to be closer to unity than to 5/3 in the optically thin 
parts of the photosphere, while in the very surface layers 
the effective Ti can become quite large (~ 8), because of 
a localized reduction of Sp. A more quantitative analysis 
of the non-adiabatic effects will require much more work, 
though, and will appear in a subsequent paper. 

Additional differences (in particular the ones reflected 
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radial mode behavior in 1-D and 3-D models, and is be- 
yond the scope of the present paper. 

However, a pre-requisite for studying mode physics ef- 
fects is a sufficient accuracy of the mean model; only if one 
includes the model effects with sufficient accuracy does 
it make sense to use remaining discrepancies to diagnose 
modal effects. 

How accurate, then, is the pressure stratification ob- 
tained from the present numerical simulations? The tests 
at various numerical resolutions show that the location, 
shape, and width of the peak of the turbulent pressure 
(relative to the total pressure) are quite insensitive to the 
numerical resolution, while the magnitude of the turbu- 
lent pressure increases slightly with increasing numerical 
resolution (Fig. ||). The magnitude scaling of the turbulent 
pressure is, however, tightly constrained by spectral line- 
broadening observations (Fig. ^|) . The existence of turbu- 
lent pressure support of about the magnitude found here 
thus cannot be doubted. Moreover, part of the elevation 
effect is due to the thermal difference between 1-D and 
average 3-D models. Any calculation of mode frequencies 
that does not include a turbulent and thermal pressure 
elevation of the upper turning points of the modes is thus 
neglecting a significant effect. If parameter fitting for such 
a calculation leads to near agreement with the observa- 
tions it merely illustrates that it is quite possible to obtain 
"the right result for the wrong reason" . 

Finally we must emphasize that while our understand- 
ing of convective effects on radial oscillations may seem 
rudimentary, it is in fact highly sophisticated by compari- 
son with our understanding of their influence on nonradial 
modes. Indeed if we consider the most nonradial mode of 
all, the f mode, in which radial and nonradial motions are 
of equal magnitude, we note that no explanation which 
seeks to replace convection modeling with a hydrostatic 
solar model can ever explain the measured frequency resid- 
uals because f-mode frequencies are largely insensitive to 
hydrostatic structure. Thus, both the f-mode frequency 
discrepancies and the remaining discrepancies in Fig. [3] 
are likely to be caused by modal effects, rather than by 
the stratification effects that we have uncovered in the 
present paper. Evidently the behaviour of both nonradial 
and radial modes needs more study. 

In order to address the modal effects on nonradial 
modes it will be necessary to invoke a more elaborate tech- 
nique than the simple planar averages we have used, a re- 
su lt already evident from th e simplified model calculations 
of Zhugzhda fc Stix (1994) . On the one hand this empha- 
sizes the gulf which still exists between theory and mea- 
surement in mode physics but, on a more positive note, it 
suggests that the remaining frequency discrepancies may 
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